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Abstract 


A new technique is presented for computing the scattering by two-dimensional structures 
of arbitrary composition. The proposed solution approach combines the usual finite ele- 
ment method with the boundary integral equation to formulate a discrete system. This 
subsequently solved via the conjugate gradient (CG) algorithm. A particular characteris- 
tic of the method is the use of rectangular boundaries to enclose the scatterer. Several of 
the resulting boundary integrals are therefore convolutions and may be evaluated via the 
fast Fourier transform (FFT)in the implementation of the CG algorithm. The solution 
approach offers the principle advantage of having O(N) memory demand and employs 
a one-dimensional FFT versus a two-dimensional FFT as required with a traditional 
implementation of the CGFFT algorithm The speed of the proposed solution method 
is compared with that of the traditional CGFFT algorithm, and results for rectangular 
bodies are given and shown to be in excellent agreement with the moment method. 
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Chapter 1 

Introduction 


Two-dimensional problems have been studied extensively from both analytical and 

numerical standpoints for many years. The demand for three-dimensional (3-D) methods 

has increased in recent years, and as a result two-dimensional (2-D) methods are finding 

their niche as testing grounds for 3-D applications. The step required to generalize a 2-D 

method to 3-D is almost always large in analytical and geometrical complexity. Most 

importantly, though, the demands in computation time and storage are often prohibitive 

for electrically large 3-D bodies. Vector and concurrent (i.e., hypercube, connection, etc) 

computers are beginning to deviate the first of these demands ([l]-[7] to mention a few 

! 

of the papers addressing this), but storage demands remain problematic. A reduction in 
storage requirements is therefore essential. 

The traditional Conjugate Gradient Fast Fourier Transform (CGFFT) method [8], 
[10] is one such frequency domain solution approach which requires O(N) storage. This 
method involves the use of FFT whose dimension equals that of the structure under con- 
sideration and therefore demands excessive computation time when used in an iterative 
algorithm. With this in mind, a new solution approach is proposed for solving scattering 
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problems that address this difficulty. The proposed method will be referred to as the 
Finite Element-Conjugate Gradient Fast Fourier Transform (FE-CGFFT) method and 
was inspired by Peters [10]. 

The FE-CGFFT method requires that the scatterer be surrounded by a double rect- 
angular box. Inside the box boundaries, the Helmholtz equation is solved via the finite 
element method. The boundary constraint is satisfied by an appropriate integral equa- 
tion, which implicitly satisfies the radiation condition. Along the parallel sides of the 
box, this integral becomes a convolution and is, therefore, amenable to evaluation via 
the FFT. The dimension of the required FFT in this method is one less than the dimen- 
sionality of the stucture leading to an 0(N ) memory demand making it attractive for 
3-D simulations. 

The proposed method is similar to the moment-method version developed by Jin 
[11]. Jin’s method was in turn based upon work published in the early 70’s by Silvester 
and Hsieh [12] and McDonald and Wexler [13], who introduced an approach to solve 
unbounded field problems. The proposed method is also similar to other methods, a 
few of which will be mentioned here. The unimoment method [14] uses finite elements 
inside a ficticious circular boundary and an eigenfunction expansion to represent the 
fields in the external region. The coefficients of the expansion are then determined 
by enforcing continuity at the finite element (FE) mesh boundary. The coupled finite- 
boundary element method [15] uses the finite element method within the boundary and 
the boundary element method to provide the additional constraint at the termination 
of the mesh. Unlike the proposed method, the solution was employed by direct matrix 
inversion as in [11], and the outer mesh boundary was not rectangular to take advantage 
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of the FFT for the evaluation of the resulting convolution integrals. Further, only one 
boundary was employed, and therefore an analytical evaluation of the Green’s function 
was required. 

In this work, we consider only rectangular structures and results derived from this 
formulation are compared to those based on traditional moment method techniques. 
Nevertheless, the proposed method is equally applicable to more complex geometries by 
using available sophisticated finite element modelling packages. 
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Chapter 2 

Analysis 


Consider a cylindrical body of arbitrary cross-section and composition illuminated 
by the plane wave 


nc (p ) = z<t>™{p) = zeikoP^B-Bo) 

(2.1) 

as indicated in fig. 2.1 (A. time dependence of e *-* has been assumed and suppressed.). 
To evaluate the fields scattered from this object, two boundaries are placed tightly around 
the body as shown in fig. 2.2. Inside the outer boundary, the Finite Element Method is 
applied to solve the Helmholtz equation given by 

V • Kp)V#/>)] + klu(p)<j>(p) = o 

(2.2) 

where 


m = E t {- P ) 

(2.3) 

v(p) = — pr 
Mp) 

u (p) = <r(p) 

(2.4) 

(2.5) 

for E- polarization and 


m = H t {p) 

(2.6) 
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Figure 2.1: Geometry of the scatterer 


v(p) = 


€ r(p) 
U (p) = Mr (p) 


( 2 . 7 ) 

( 2 . 8 ) 


for H-polarization. Also, k 0 = is the wave number, and p r and t r are the relative 

permeablility and permittivity, respectively. 

The appropriate boundary condition is enforced on the surface of the impenetrable 

boundary, while the radiation condition is satisfied implicitly by evaluating the integral 
equation 

m = - £ {<%?> [£*«] . *on } a (2.9, 

on the boundary T a , where 


<?(?.?) = 


( 2 . 10 ) 
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Figure 2.2: Partially discretized body 

is the 2-D Green's function in which *<’>(.) denote, the zeroth order Hankel function 

of the second kind. Furthermore, f and Jt are the usual observation and source position 
vectors, respectively, and 

\P-t\ = \/(* -*')’ + Of -yO 1 (2.11) 

The normal derivative, are taken in the direction of the outward normal of r c . 

2.X Discretization of the Object and Field Quantities 

In fig. 2.2, r. is the field/observation point boundary, and T c is the integration 
contour, which is placed midway between T. and f,. Also, 1. denote, the contour 


9 




Definitions for Finite Element Mesh 
N g = total number of nodes in the finite element mesh 
N e = total number of elements in the finite element mesh 
N x = number of nodes on T a or along the i-direction 
N y = number of nodes on T 0 or T*, along the y-direction 
N a = total number of nodes on r o 
Nb = total number of nodes on Tj 
Nab = N a + N b 

r a = ZUi r a , 
r k = E? =1 T bi 
r, = EU t Cj 

( x a, , Va , ) - coordinates of a point on contour r ai 
( x b,i y b , ) - coordinates of a point on contour T bi 

Table 2.1: Defimtions for the finite element mesh 

enclosing the impenetrable portion of the scattered Herewith, each side of T a , T b and 
T c will be numbered counterclockwise starting from the top side, as indicated in fig. 2.2. 
The fields in the region between T a and T d satisfy (2.2) in conjunction with the required 
boundary condition on T d . The boundary integral equation (2.9) will be enforced on T a . 

To numerically solve (2.2), it is required that the region within T a be discretized. 
This is done in a traditional manner when employing the finite element method. The 
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Definitions of Field Vectors (in terms of field unknowns at nodal points) 


4> a , - fields corresponding to the nodes on r o . 

K = MdS corre8 P on ‘^ n 8 to the first N . o, N, (whichever is appropriate) nodes on T ( , 
4>b, = fields corresponding to the nodes on 

4>i = fields corresponding to region I enclosed by r 6 and 
<fa = fields corresponding to the nodes on the T* 


Table 2.2: Definitions of the field vectors 

globai node numbering starts from the right endpoint of contour r., and proceed, coun- 
terclockwise. The numbering continues beginning at the right endpoint of contour T,, 
and proceeds counterclockwise. Within T,. the node, are numbered from the lower teft 
corner and proceed column by column. The definition, pertaining to the FE mesh are 
given in table 2.1. Each node corresponds to an unknown field value to be determined. 
It is important to associate the unknown field value, corresponding to the various node 
group, on contours T. and f, hy using different variables The labeling scheme is given 

in table 2.2, and this discrimination of the nodal field, is required, since they are treated 
differently in the analysis. 


2.2 Derivation of the Finite Element Matrix 

One of several approve, may be used to generate the finite element matrix, such a, 
the variational approach or the method of weighted mridujs. In this development, we 
wiU utilise Galer kin’s method, which is a special ewe of the method of weighted residuals 
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with the distinction that the testing functions are the same as the basis functions. 


Proceeding with the finite element analysis, we may rewrite (2.2) as 


dx 


v(x,y) 


d_ 

dx 


<t>(z,y /)] 


+ 


d_ 

dy 


v(x,y)—<t>(x,y) 


+ k 2 0 u{x , y)<j>(x, y) = 0 (2.12) 


the residual of which is given by 
' 9 


d d * d r ^ 1 

dx v(<X,y ^dx^ X,y ^ ~ dy [ t '( a! ’ y )^^ x ’ y )J “ (2.13) 

The field within T a may be represented as a summation of piecewise continuous functions 
and, thus, may be written as 


1V« 


d>(x,y) = £<£'(*, y) 

e=l 

where <j> e (x, y) is the field within element e. It is expanded as 

<£'(*, y) ~ £*/(«, F #s 


(2.14) 


j=i 


(2.15) 


where N ‘\ s are the standard shape functions (found in any standard finite element book), 
’s are the fields at the nodes, and n is the number of nodes per element. The weighted 
residual equation for the eth element is defined by 


Jl NfRdxdy = 0 t = (2.16) 

where S e denotes the surface area of the eth element. Inserting (2.15) and (2.13) into 
(2.16) yields 


tjL 4iHSK(-f)- 


kluN* 


j 4> e jdxdy = 0 
i — 1,2, ..., n 


(2.17) 
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Further, by using the identities 



dNf dm 
dx dx 


dNf dm 

v 5 

dy dy 


and the divergence theorem 


(2.18) 


(2.19) 


JL [tx + 1 ;) dxiy = i. + <*> ■ Ml 

where C e is the contour enclosing the eth element, (2.17) becomes 

A / r ( dm dm dm dm \ 

E// s . + sf-sf - W* 

= E f c _ K ^ (.-gA. + v-^-y j.ndl i = 1, 2, .... n 


( 2 . 20 ) 


( 2 . 21 ) 


This may be written in matrix form as 


A e <f> e = b e 


( 2 . 22 ) 


where 


and 



(2.23) 


(2.24) 


For linear triangular elements, jV* are given by 






and 


fi e = 


1 , 
-del 


i y\ 

1 x\ y| 
1 *3 2/3 


= -(bfc' - btf) 


(2.26) 


a. = 


6 ? = 


c- = 


z'j/Jt - 4 j/> 


y' - vl 

xl-x‘ 


(2.27) 

(2.28) 
(2.29) 


where (z*, y*) is the coordinate of the ith node of the eth element. From (2.25) 


dNf bj 
2 fi e 
c? 

2fi« 


0z 

diVf 

01/ 


(2.30) 

(2.31) 


Substituting these and the formula 


Us 


s .wmi 


into (2.23), we find 


mt# - + c?cj) - *;«•£«« 


(2.32) 


(2.33) 


where 


= 


2 if t = j 
1 otherwise 

In (2.33) we have assumed that u and v (the material constitutive parameters) 
constant within each element and are given by ti e and v e , respectively. By summing 


(2.34) 


are 


over 
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an elements as implied by (2.14), we may write the overall system more explicitly as 


■r4<ja A a fr 0 0 


4>a 


b o 

Aba A^ Abi 0 


<t>b 


0 

0 Aib An Ai& 


4>i 


0 

0 0 An A dd 


<i>d 


0 


In this, the values of the elements in the submatrix A pq are the contributions associated 
with the nodes in group p which are connected directly to the nodes in group q. 

One can easily show that the line integral contribution (see 2.24) of those elements 
vanishes everywhere, unless the element has a side in common with T a . As a result, 6' 
contributes only from the boundary T a of the finite element region, as indicated by the 
presence of the vector 6 a in (2.35). Without a priori knowledge of the total field on that 
boundary, 6“ cannot be determine. We may, however, provide the appropriate condition 
on this boundary by utilizing the integral equation (2.9). This amounts to replacing 

the first block-row of the matrix (associated with <t> a on T a ) with a discrete form of this 
integral equation. 

2.3 Evaluation of Boundary Integral 


The boundary integral in (2.9) may be written ae a enmmation of four integrals, one 
for each side of the contour T c as 

m = *"“(« - [<?»,?') -g^m - m-gjG®,?) J 

+ /r„ + 


15 






where the derivatives along the x and y directions, denoted by 5 ^ and 3 ^, respectively, 
have been left in this form for the later convenience of determining them numerically. 
More explicitly, we have 


4>(x,y ai ) = 4> xnc {x,y ai ) 


- | J r C?(*-x',y a ,,y Cl )~^(i / ,y cl )-^(i , ,y cl )-^ ; G(x-i / ,y ai ,y cl )| dx' 


dn y ’* CI/ irv ~’* c l, dy> 


+ j r \g(x - x \y ai ^y Ci )^4>{x\y Ci ) + (t^x^yc^-j^Gix - x ',y ai ^ ,y C3 )j dx' 

+ j r ,y / )^<K*c4,y') - <A(®c 4 ,y')^ 7 G ; (*»a;c 4 ,yai ,y')] <*y'} 


(2.37) 


and 


4>{x ai ,y) = <t> inc (x a 2 ,y) 


0 3 

G(x a2 ,x,y, j/cj ) 4 >(x , ye, ) ~ <A(x , y Cl ) ^6(1#, 1 1 1 J/i ye, ) 


dx' 

<V 


-(x. 

/ r ^ ^ 

> ^ca ? y ~ y > 1/) + > y ^ y “■ y ) 

* r 3 3 

+ j r iGixa^x'.y^^—^x^y^ + ^x^y^—Gixa^x^y^c:,) 

+ J r [<?(*aj,x Ct ,y-y / )^-^(x Ct ,y / )-^(*c4,y , )^;G(xa,,Xe 4 ,y-y') dy'j 


(2.38) 


where the first subscript on x or y refers to the contour (a, b or c), and the second refers 
to the contour number (see fig. 2.2 and table 2.1). It is noted that the arguments of the 





Green’s functions have been modified to imply a convolution when appropriate, and this 
representation will be used throughout. 

The normal derivatives of <f> may be evaluated via the central difference formulas 

Q^^( x c,y') = — [4>{xa,y') - <t>{x b ,y')\ + 0(A 2 ) (2.39) 

d , i 

,yc) = — [<t>(x',y a ) -<Hx\y h )) + 0(A 2 ) (2.40) 

where A is the displacement of T a from T* (A is usuaUy less than one tenth of a wave- 
length). Substituting (2.39) and (2.40) into (2.37) and (2.38), we obtain 


<t>(x,y ai ) = <f> inc (x,y ai ) 

3 3 

" {l ei [ K;G{X ~ fcl ) - K V + G(* - f, J/ai , y Cl )<Kx\ ft, )] dx' 

+ Jr C3 \ K Z G ( X ' x <* » 1 . y')<K*a % , y') - K~ G(x, x C3 ,y a3 , y')<j>(x h , y')j dy' 

Ir C3 \^ v G ^ X X ’ ^° i ’ 3 ) < t > i x » y» 3 ) — K y G(x — x', y„ j , y C3 )<j>{x ' , j /j 3 )j dx' 

+ Ir Ci G ( x ' Xe *'y*i >y ) < K x niy ) “ K^G(x,x C4 ,y a3 , j/VOr&i , j/)J dy'| 


(2.41) 


and 


<H*a 2 ,y) = <t> inc (x ai ,y) 

4 4 

{lr ci [ ir v <sr ( x «j.*'»y.»c 1 )^(* / ,y« 1 )-^+C?(x a ,,* , ,y,y cl )0( x ' ty6i )j < i a . / 

+ I rct [** G (*“2 . *c, , y - , y') - K~G(x a3 ,x Ci , y - y')<j>(x h , y')J dy' 

+ L Ci \ K v G ( Xa \ * **» *» )#*', y as )~ K~ G(x a , , x', y, Vc3 )(^x\ y bi )j dx' 

+ Ir u \ K * G ( X *\ ' Xc *>y- , y') - K+ G(x., , x C4 , y - y')<t>{x hi , y')J dy'J 

(2.42) 
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in which 


*?« 

K v = 


i ± IA 

A 2dx' 

1±IA 

A 2 By’ 


(2.43) 

(2.44) 


Assuming a pulse basis expansion for the nodal fields (centered at the nodal positions 
along contours T a and 1^), a midpoint integration may be performed for the evaluation 
of the integrals in (2.41) and (2.42), to obtain 


) = <t> ,nc (xi,y ai ) 


[ K v x ^yai,yci)4>(xj,yai)- K+G{xi - Xj,y a i,y Cl )<t>(xj,y bl 

^ ^2 i*' G(Xi,x e3 ,y at ,yj)<t>(x ai ,yj) - K~G(xi,x Ci ,y ai ,y, )</>(£&,, y,) 

j = 1 1-3 3 

\ KyG(x t - ij,y ai ,y C3 )<t>{xj,ya i ) - K y G(ii — x j,y ai ,yc 3 )4>(xj,yb 3 ) 

j=l L 3 3 

+ A ^ G(ii , x Ci , y„ , , yj)<t>(xa4 ,yj) - K+ G(xi, x ^ , y„ , , y> )<£(x6 4 , yy ) > 

>=i L 3 3 J I 


(2.45) 


and 


<K x a,,Vi) = <£ mc (*«a,y») 


l ^ Xrf k >*i> Vi^Vci )<K x i,y ai ) — K y G(ia 3 , *j,y,, y Cl )4>{xj,y hl 

( i- 1 1 4 4 

■*" ^ ^( x « j i x ca * yi - y> )^(* aa * y> ) - AT~ <?(i a , , Xc, , y< - y, )<t>(xb 3 , y, )j 

+ ^ X* G{x a ^ ,Xj,yi,y C3 )4>(xj,y ai ) — K y G(x aj ,Xj,yi,y Ci )(p[xj,yb 3 )^ 
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r 


+ A 12 [ K X G(x a , , *e« , y, - y : )4>(x at , y, ) - Ff + G(x a 2 , x C4 , y, - Vj )4>(x bi , y } ) 
i= 1 14 4 


(2.46) 


In these x, and y, denote the ith matching/testing points corresponding to the nodal 
locations on T a , while Xj and yj denote locations on IV We recognize some of the terms 
in (2.45) and (2.46) as discrete convolutions amenable to numerical evaluation via FFT. 
The subsystem (2.45) and (2.46) may be written more concisely as 




K c 


S a,b! 

T U 

Kb 3 R * 

1 

if 


K 

4 > d2 


K c 


Kb i 

Kb, 

Kb 3 

^a,b. Ry 



^03 


6 inc 

V-'aa 


^o 3 6 , R * 

K* 

S a 3 b 3 

Kb. 


K 

<^4 


4> inc 

rat 


Kb , 

Kb,*v 

Kb 3 

Kb. 


$ a « _ 





Kh 

S aib* Rx 

Kb. 


K 



+ 



T~ 

aibz 

Kb.Ry 


K 




S UB. 

Kb 

Kh 

K^ 


K 




.K>> 

S •**»&* 

Kh 

£ + 
i 


. K . 


(2-47) 

with the various parameters to be given explicitly later. The matrices R XiV simply re- 
verse the order of the unknown vector so that the convolutions may be performed prop- 
erly. This is required solely because of the employed counterclockwise nodal numbering 
scheme. 

Since 

Wi« , Wto t f i * l ’ U4 <**> 

J = 2,3,4,! 
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the vector 


$>i fib-t K ^64 ] 


(2.49) 


can be related to the actual unknowns on the contour Tj via a transformation as 


4>'b = Db4>b 

and (2.47) may then be written as 

(/ + L aa )4> a - L a bDb4>b = 


(2.50) 


(2.51) 


or 


where 


and 


1 4* L aa — 


L a b — 


1 Bl 


<P a 
4>b 


= K nc 



K* 

S tb, R * Kb. 


T~ 

*2*1 

l + s ih T U 

*^0364 Ry 


Kb, 

1 + Kb. 


T at bi 

st*** 


! + Kb t \ 

s i>, 

Kb> 

Kb,** 



<r+ 

*2*1 


^*363 



s u*' 

T z* 

Kb, 

Kk 




Kb, 

S U 



(2.52) 


(2.53) 


(2.54) 
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After replacing the first block row of (2.35) with (2.51), the complete system may be 
finally written as 


f "h T ao L ab D b 0 0 


4>a 


K nc 

Aba A bb A(,/ 0 


4>b 


0 

0 A/j An A/j 


<t>i 


0 

0 0 Adi Add 


4>d _ 


0 


to be solved via the CG algorithm. 

The Ele ments of A fl/ 

The elements of A BI defined above may be evaluated via the discrete Fourier trans- 
form. Specifically, we have 

^ 6, = DFT 1 [DFT[G{x,y ai ,y hi ) ± G y (x,y ai ,y bl )]DFT[d> tl ]1 (2.56) 

3 ^ 3 3 3 J 

= DFT ~ 1 { £>FT [ G ( I »l'«3*yfci)±^»(a:,ya„y6 1 )pFT[^ I ]l (2.57) 

3 k 3 3 3 J 

= DFT 1 (DFT[G(x ai ,x b2 ,y) ±G x (x at ,x hl ,y)\DFT[<t> hi )\ (2.58) 

4 V 4 4 4 J 

= DFT~ l |£)/T[G(z a4 ,x 6 , ,y)±G x (x ai ,x b ,,y)]DFT[<t> b ,]^ (2.59) 

in which DFT denotes the discrete Fourier transform operator. Also 

G{x,x\y,y') = ^H^(k 0 \p - p'|) = I±H™ (k 0 J(x - x'f + (j, - j,') 2 ) (2.60) 
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G x{z,x' ,y,y f ) _ — — Q^x^fyyt^ 

8 v'Cz - x') 2 + (y - y') 2 

c,(*, *',»,»') = j-^pG(x,x\y,y') 


(x - X ') (2.61) 


8 v /(z-x') 2 + (y-y') 5 vy 3 

Special cases of the convolution operators for the chosen mesh are given as 
n . ; j ^ (k 0 J&Z — xj 3 )** + (y — y') 2 ) 

_ .It. . w . w ; ^ — 4-( — \ A U V 4 4 f \ 


(y - y') (2.62) 


Gx(xa 3 ,x bi ,y-y') = ±(-±)Ak 0 - 

<4 O 


(^°2 - s * 2 r + (y - y 7 ) 

4 4 


1^02 X b 7 

4 4 


{*;; ( 2.63) 


r , , s > ~ 1 ) +(»«i- y>, ) 2 ) 

G v (z-x,y ai ,y 6l ) = ±(_£)A* o - v = ' J » 7 In _ Vfc I 

3 3 8 J(* -*')* + (»«, -y *,) 2 |y °i y6 r 

V 3 3 

and the corresponding expressions for G are implied by the arguments in the previous 
two equations. 


The cross-term element submatrices are given by 


K«i»aj. ~ G{Xi,x b% ,y*, ,y,) ± <?,(*,, XJ 3 ,y„, ,yy) (2.65) 

L 3 4 j I; 43 43 

l^aa&i I .. = ^(*«a ^jiViiVbx ) ± G v (x a3 >z>.yi»y»i ) (2.66) 

1 « 3 J V 4 3 4 3 


Gx{%\y 2b 2 y V*\ yjfj) 


_ ±(Z i )k g l 2 ) (*o^/( x » - ) 2 7 (yq, - y>) 2 ) 


V 1 y \ * 3 ’ ) ( Xk. 

- “ ) J + (*., - y,) ! , * i_ **S |A {*,. 
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and 



where again the corresponding expressions for G are implied by the arguments of the 


previous two equations. Making the substitutions 


and 


(x,-x fc ,) 2 = (*-.5) 2 A 2 (2.69) 

4 

(V.. - K? = 3 ’A’ (2.70) 

3 

= A*/(i-.5 yT? (2.71) 



(x OJ ~ Xj) 2 = j 2 A 2 (2.72) 

4 

(*i -»,)’ = (i-.5) I A 3 (2.73) 

3 

=>y(x a j -x J ) 2 + (y,-yfc,) 2 = A^i 2 + (*-.5) 2 (2.74) 

we may write G r and <7 y as 



to be used in the actual implementation of the system. Since each of the above relations 
are similar, we are required to store only one of them and alter the signs accordingly. It 
should be noted that, however, this is not the most efficient method of storage. Storing 
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only a few of the cross term values and using an interpolation scheme will reduce the 
storage considerably. Of course, an interpolation table of (2.75) and (2.76) will lead to 
a substantial reduction in memory at the expense of some computational efficiency. 
Assuming that the positive sign is chosen in equations (2.75) and (2.76), we have 


T+ — q — 

a 2 fcj a 2 bi 

4 4 

T+ — T+ 
a i h *ih 

3 3 

T+ — 'T+ 

0 2^3 ^0363 

4 4 

T+ — 7— 

a 1 64 a 1 64 

3 3 


(2.77) 


Choosing the positive sign for the (2.63) and (2.64), we also find 


3 3 

^ b 2 ^<*2 b 2 

4 4 

^*3*1 ~ ^a 3 6i 

3 3 

S t A b 2 = 2 

4 4 


(2.78) 


Thus, the elements of A Bl may be written as 




St.*,** 


% 

1 + s :,* 

% 




< + 





’ + 
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Lab — 



T a x b? 

S ai b 3 Rz 

^164 

7 — 

a 2 &i 


7 — 

0363 

^a 3 b t 



S U, 

■^364 

■^4 <>1 

h R v 




( 2 . 79 ) 


The elements of the adjoint of A SI required in the implementation of the CG algorithm 


are 


(I + L aa ) a = 


( L ah D h f = 


DT 


r + («.-,». r 

at*,)* 


at*,)* 

(T* h Y 


at*,r 

*fa + .*,)“ 


at H r 

'+at*,)“ 

at»,>“ 



at..r 

'+a-.*.)“ 

at*,)* 

at*,)” 


at*,)” 

at*,)' 

at*,r 

at*,)” 

*f(W 

*1(5:,.,)' 

' at*,r 

(■£*,)“ 

at*,r 



at*.) 1 

(•s;*.r 


(2.80) 


2.4 A CGFFT Algorithm 

The CG algorithm to be employed for solving the system (2.55) is as follows 
Initialize the residual and search vectors 


7fc = II [ 0 0 0 ] T IlHl II 2 

s = A<t>(°) 


(2.81) 

(2.82) 
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6-3 


= 6-3 

a = A a rW 

(2.83) 

II 

(2.84) 

/* 01 = 7,- 1 

(2.85) 

pf) = 

( 2 . 86 ) 

Iterate for fc = 1,..., jV 3 

5 = ApW 

(2.87) 

7. = |M|? 

(2.88) 

a (/t) = 77 1 

(2.89) 

*(*+*> = /) + a <y) 

(2.90) 

r (i+l) _ r (*) _ a (fc)p(*) 

(2.91) 

7r = II II’ 

(2.92) 

3 = y4“r (<t+1) 

(2.93) 

7# = || * 111 

(2.94) 

7 « 

II 

jT 

(2.95) 

p(*+!) = p(*) +/ 3<*) 3 (*) 

(2.96) 

Terminate when k = N g or < tolerance. 

The individual operations associated with the A BI matrix-vector products are quite 

numerous and, therefore, will not be shown explicitly. However, it 

can be shown that 


the total system may be decomposed into a summation of two matrices; one involving 
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operators associated with the boundary integral and another involving the elements of 
the finite element matrix. The system matrix may then be written as 


W — + {sfe}- 


where 


and 



I + L 

aa 

LabDb 

0 0 

2\ 

{*B/} = 

0 


0 

0 0 

22 


0 


0 

0 0 

23 


0 


0 

0 0 

24 


0 

0 

0 

0 

2\ 

{^fe} = 

Aba 

Abb 

Abi 

0 

22 


0 

Aib 

An 

Au 

23 


0 

0 

Adi 

Add 

24 


and 


{**/} = 


{*fe} = 


" 


- 

o o o 


2\ 

DlL\ b 0 0 0 


22 

0 0 0 0 


23 

0 0 0 0 


24 

r 



O 

r* 

o 

o 


2\ 

o 

S* 

Ju 

o> 

O 


z 2 

0 A bI AJj A%j 


23 

. 0 0 A% _ 


24 

L J 


( 2 . 97 ) 


( 2 . 98 ) 


( 2 . 99 ) 


( 2 . 100 ) 


( 2 . 101 ) 
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In each case, the operation is performed such that only the resulting vector {s} need 


be stored. 
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Chapter 3 

Computational Considerations 


This method is efficient in terms of memory usage and computation time. We discuss 
each of these aspects in detail below. 

3.1 Storage Efficiency 

The fundamental advantage of this method is the reduction of storage requirements, 
so that the scattering by electrically large bodies may be evaluated. To show that the low 
storage requirement of 0(N g ) is assured, we refer to tables 3.1 and 3.2. These contain a 
list of aU major memory consuming variables. A summarized list is also given in table 3.3. 
Specifically, table 3.3 includes the memory requirements pertaining to the finite element 
matrix (FE), fast Fourier transforms (FT), boundary integral cross terms (Cross) and 
conjugate gradient variables (CG). N c is one less than the number of elements connected 
to a particular node, and a typical value of 5 is used here. 

To put the quantities of table 3.3 in terms of N„ the total number of nodes, we 
consider two special geometries. The mesh in fig. 3.1 corresponds to a penetrable 
body, while that of fig. 3.2 corresponds to an impenetrable structure tightly enclosed 
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by the picture frame. Within each special case two extremes are considered; a mesh 
corresponding to a square object (N x = N y ) and a long strip ( N x >> N y ). In each case, 
N g is assumed to be large. 

Alluding to table 3.4 the total storage is 0(N g ) for the square region, but is some- 
where between 0(N g ) and O(Ng) for the ( N x >> N y ) case. This is based on the 
assumption that all cross terms are individually stored, but by using an interpolation 
table, the 0(N g ) memory requirement can be assured regardless of the value of N x with 
respect to N y . In table 3.5, more dramatic results for the storage of the cross term are 
listed. In this case, all of the unknowns are on the outer two boundaries, so it appears 
that the storage is 0(N *) for the square case. One must note, however, that the factor 
multiplying the N g term may be quite small. The strip case, on the other hand, requires 
an O(Ng) storage. This case would be am unlikely candidate for the use of this method, 
since that structure would be handled much more efficiently via a direct implementation 
of the CGFFT method. As noted above, the storage of the cross terms may be brought 
down to 0(N g ) for all cases by using am interpolation table, and this will certainly be 
necessary in a 3-D implementation. 
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Memory Consumption 


variable 

type 

count 

comment 

Mesh 


Xg 

R 

Ng 

X coordinate of global nodes 

Yg 

R 

*9 

Y coordinate of global nodes 

Nglob 

I 

3N e 

Needed for FE matrix formation 

Pointers 


ABint 

m 

N ab 

Points to observation and integration points 

Pnodes 

■ 

Pnum 

Points to nodes on conducting bodies 

dmatl 

■ 

N g - N ab 

Element material properties 

Finite El< 

sment 

Matrix (FE) 

Ar 


~ i^XNg - N a ) 

Non-zero values of FE matrix 

col 

i 

~ - N a ) 

Column pointer of FE matrix 

[ row 

i 

N g -N a 

Pointer to rows of FE matrix 


Conjugate Gradient (CG) 


Phiz 

C 

N g 

Unknown vector 

CGI 

C 

Ng 

Residual vector 

CG2 

c 

Ng 

Search vector 

CG3 

c 

Ng 

Temporary vector 

q 

c 

MAX(N t ,N y ) 

Temporary vector 

phiinc 

c 

N a 

Incident field vector 


Table 3.1: List of major memory-consuming variables 
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~ UAX(N z ,N y ) 


Legend 
R = REAL*4 
C = COMPLEX 
I = INTEGERS 


Table 3.2: List of major memory-consuming variables (continued) 






Figure 3.1: Example of the mesh of a penetrable structure 



Figure 3.2: Example of the mesh of an impenetrable structure 
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Major Memory Consumption (N x > N y ) 

Item 

Type 

Count 

FE 

COMPLEX 

(£L±l )[Ng _ 2 (N x + N y )\ 

FT 

COMPLEX 

12 N x + 8 N y 

Cross 

COMPLEX 

2 Nl 

CG 

COMPLEX 

4Ng 


Table 3.3: Summary of major memory consumption 


Major Memory Consumption: Penetrable 

Item 

£ 

II 

N x » Ny 

FE 


( £L f 1 )Ng( 1 ~ jfe) 

FT 

20^ 

12Ng/(Ny + 2) 

Cross 

2Ng 

2 (Ng/Ny)> 

CG 

4Ng 

4Ng 

total 

~ 9Ng 

~ 2 (l^+ j) 2 + 6 2T1^ + 7N 9 


Table 3.4: Summary of major memory consumption: filled mesh 
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Major Memory Consumption: Impenetrable || 

Item 

ii 

N X » Ny 

FE 

(&±l)N g /2 

(&±l)N g /2 

FFT 

5N g /2 

3N g 

Cross 

A^/32 

NJ/S 

| CG 

4 N, 

4 N, 

total 

~ N 2 / 32 + 8 N g 

~ + llNg/2 


Table 3.5: Summary of major memory consumption: open mesh 

3.2 Computational Efficiency 

Since the primary importance of this method is storage reduction, a comparable level 
of efficiency with alternative methods is a bonus. A method for which a fair comparison 
may be made is the standard CGFFT. This requires a 2-D FFT, which is slower than 
using multiple 1-D FFTs for large bodies. We compared the two methods for a specific 
penetrable scatterer using an , -Apollo 3500 without code optimization. The pertinent 
CPU times are compared in table 3.6. The comparison provides only a relative measure 
of the speed difference. 
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Body Properties 

FE-CGFFT 

CGFFT ( 

Composition 

Dimensions 

T/I (s) 

I 

Total 

T/I (.) 

I 

Total 

dielectric 
(r = 4 - j. 1 

2A x 2A 

8.6 

155 

1333 

170 

33 

5610 


Legend 


T/I = time/iteration 
I = number of iterations 


Table 3.6: A comparision of computation efficiency of the FE-CGFFT with the CGFFT 
method 
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Chapter 4 

Far Field Computation 


The scattered fields may be computed as 

= ■ i. \t'^\ - [tM-p,?)) } # (4.1) 

Using the discretization scheme developed earlier, we have 

4> 3 (x,y) = 

{i ei \ K V G ^ *'> Vex )<K*\ ) - K+G(x, x', y, y Cl )#*', y bl )j dx' 

+ Jr Cj ^ K * G ( x,Xc i ’ y* y')<K x *2 » lO - G(x, ic, , y, , y')] dy' 

+ J rej [■Ky'G(x, x' , y, y^ y a3 ) — K~G(x,x r , y, y^ )<j>(x\ y bj )j dx' 

+ Jr c (^r <?(*,**, y.y'Mxa^y')- K+G(x,Xc4,y,y')<t>(x h ^y')} dy'\ 


where the definitions for K * and A"* 


: are 

as specified previously. Letting 

(4.2) 

1 , 

/•*>+$ „ . 


A J, 

G{x,x',y,y c )dx' 

Z J~2 

(4.3) 

1 i 



A l 

G(x,x c ,y,y )dy' 

»j-T 

(4.4) 
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(4.5) 


1 f x J+ 2 d 

7 (x,y,y c ) = 2 J x _a lh/ G ( x,x '' y ' yc ' >dx ' 

1 d 

7 y ( x,x c ,y ) = - y _ a — G(*,x c ,y,y')dy' 


(4.6) 


(4.2) becomes 


4> a {x,y) = 

I Nx 

- < E (&**(*. ) - 7 X (x,y, Vex )] {K h - [0 x (x, y, Vex ) + 1 X (X, y, y cx )] {&, } } ) 


li=i 

Ny 

+ + 7 1 '( z > a; c 2 ,2/)] {<^o 2 }> - [(3 v (x,x C2 , y) - 7 y (x, x C2 , y)] {^{, 2 }j) 

j= 1 

JVx 

+ E ) + 7*(*, y, y C3 )] {4> a3 }j - [/? x (x, y, y C3 ) - 7 x (x, y, y C3 )] {0 6j }_, ) 

j=i 

Ny 

+ E (10* (*» x « ,V)~ l y (x, x Ci , y)] {4> a< }j - [/? v (x, x Ct , y ) + 7 ®(x, x C4 , y)] {<& 4 },) 


j=i 


(4.7) 


we 


valid for all observation points (x,y). The specialize (4.7) to far zone computations, 
must introduce the appropriate asymptotic expansion for the Hankel functions implied 
in (4.3)-(4.6). In doing so, we have 


where 


0 x (x,y,y Ci ) = j7o(/>)/i(d,y cl )^ co - fl 

3 3 

P v (x,x C3 ,y) = jf o (p)f 2 (0,x cj )e^“«* 

* 4 

7 X {x,y,y ct ) = -f o (p)fi(0,yc 1 )k o Asin0e jk ° x ’ cote 

3 3 

7 v (x,x cj) y) = -f o (p)M0,x Cl ) k 0 Acos$e jk °'>’ tin6 




(4.8) 

(4.9) 

(4.10) 

(4.11) 


(4.12) 
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(4.13) 


jfcoyc, sin 9 

/x(^.y Cl ) = -e 3 

3 

j*0*C 2 oosfl 

h\"lXc2 ) = -e 4 

4 


(4.14) 


in which (p, 0) imply the usual cylindrical coordinates of the observation point. Substr- 
ing expressions (4.8)-(4.11) into (4.7), we obtain 


x,y) = -Up) 

( 

(U + *oA Sin 6] { <p ai -\j- k Q A sin 0) {4> bl } j )f l (0 , y Cl )e jk ° x > co * e 


U=l 

+ ]C CC7 - *oA cos 0} {<f> aj }j - [ 7 + cos 0] }_,) / 2 (0, sin * 

i=i 
N x 

+ " ^Asin 0} {U}j -\j + * o Asin0] {^ 3 } J )/i(<?,j/ C3 )e^ C089 

i=i 
N v 

+ H (L/ + *oA cos 0] {<^ 04 }j-U - k 0 A cos 0] {& 4 }_,) / 2 (0, I C4 )eJ*»w ‘in « 
i - 1 


The echowidth is defined by 


and from (4.15) we have 


\4>U 

a = lim 2 ttp , . 

« * inn 


P— oo “ |<^ ,nc j2 


(4.15) 


(4.16) 


a = 


4£ 0 


N, 


£([7 + *oAsin0]{^ ai }_, - [j - k 0 A sin 0} {&,},•) f x {0, y Cl )e> k ° x > co * 6 
i=i 

N y 

+ S(D'“ *<>Acos0] {<^ aa }j - \j + * 0 Acose]{^ 62 } > )/ 2 (^,x C2 )e J ^ 8intf 

i=i 
Nm 

+ 2(b'- ^A sin 0] {0 a3 }j ~[j + k a A sin 0] {<&,},) ^(0, y C3 )e jfc ° x ^ C08 * 

i= i 

% 

+ £ (U + *°A cos 0] {<fc, 4 } ; -[j- k 0 A cos 0] {<t> h }_,) / 2 (0, Xct )e jk °'>> “ 5 

j=i 


(4.17) 
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Chapter 5 

Code Validation 


The scattering patterns from several rectangular structures are presented. The echo 
width is computed for each structure and compared to the results of the moment method. 
The bodies are discretized at a sampling rate of 20 samples/free-space wavelength. 
Results are presented for the following cases: 

• perfectly conducting bodies (figs. 5.1 and 5.2) 

• partially and fully coated perfectly conducting cylinders (figs. 5.3 - 5.8) 

• composite body (fig. 5.9) 

In each figure, the discretized geometry is shown along with the corresponding results. 
As seen in all cases, the generated patterns using the FE-CGFFT formulation are in 
excellent agreement with the moment method data. 
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FE-GCFFT 


MOM 







Angle (deg) 


Figure 5.3: E z backscatter from a .25 x 1A perfect conductor with a A/20 thick material 
coating containing the properties e r = 5. - j.5,fi r = 1. 
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.25 x 1. X Coated Perfect Conductor (H-pol) 



Angle (deg) 


Figure 5.4: H z backscatter from a .25 x 1A perfect conductor with a A/20 thick material 
coating containing the properties e r = 5. - j. 5,/i r = 1. 







Angle (deg) 


Figure 5.5: E z backscatter from a .25 x 1A perfect conductor with a A/20 thick material 
coating containing the properties c r = 5. — j.5,/i r = 1.5 — j. 5. 
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.25 x 1 . X Coated Perfect Conductor (H-pol) 



Figure 5.6: H z backscatter from a .25 x 1A perfect conductor with a A/20 thick materia] 
coating containing the properties e r = 5. - j.5,/i r = 1.5 - j.5. 
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.25 x 1. X Double Topped Conductor (E-pol) 



Angle (deg) 


Figure 5.7: E z backscatter from a .25 x 1A perfect conductor with two A/20 thick top 

material coatings. The lower layer has the properties e r = 2. - j.5,/x r = 1.5 — j. 2, and 

the upper layer has properties e r = 4. - j. 5,/i r = 1.5 — j. 2. 

47 







.25 x 1. A. Double Topped Conductor (H-pol) 



Angle (deg) 

Figure 5.8. H z backscatter from a .25 x 1A perfect conductor with two A/20 thick top 
material coatings. The lower layer has the properties c r = 2. - >.5,/i r = 1.5 - j. 2, and 
the upper layer has properties e r = 4. - j.5,n r = 1.5 - j. 2 . 





Backscatter Echo Width/* (dB) 


.5 x 1. X Composite Body (H-pol) 



Angle (deg) 


Figure 5.9: H z scattering from a composite body. Both the perfect conductor and the 
dielectric body are A/2 in each dimension. The material properties are e T = 5.-j.5,^t r = 







Chapter 6 

Conclusions nnd Future ^Vork 

A procedure was developed for computing the scattering by 2-D structures. This 
procedure combined the boundary integral and finite element methods, and the result- 
ing system was solved via CGFFT. The main advantage of the new approach was a 
reduction in memory demand to O(JV) compared to 0(AT>) required with traditional 
solution techniques. A detailed map of the storage requirements was presented, and the 
principle memory consuming arrays were discussed. Also, the computational efficiency 
of the technique was examined and found favorable. To validate the proposed solution 

approach, several backscatter patterns were presented and compared with results based 
on traditional solution methods. 

A goal is to extend the technique to 3-D applications. In this case, the cross terms 
must be efficiently stored using an interpolation table to ensure an O(N) storage re- 
quirement. Also, the use of a simple boundary (as in [15]) in the application of the 
boundary integral equation would be desirable for additional storage reduction. Higher 
order elements are further of interest to increase the CG convergence rate. Second order 
elements are also within the solution domain of Maxwell’s equations and would allow 
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a more accurate evaluation of the normal field derivatives. In addition, there are other 
numerical difficulties that must be addressed in 3-D applications. The modeling of the 
fields near corners of the scattererrequires some care (an obvious approach is to avoid 
placing a node at the corner location). Also, the field discontinuity at material transi- 
tions must be handled properly. The standard basis functions ensures continuity across 
a boundary, but this will require some modification. 
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Chapter 7 

Program Manual 


In this section, brief descriptions of the pre-processing programs (mesh generators) 
and the main processing program are provided. They have been executed on an Apollo 
workstation an IBM 3090-600E and a Cray Y/MP. 

7.1 Description of FECGFFT 

The main processing program, FECGFFT, is a menu driven program which allows 
the user to load the desired pre-generated mesh file, choose the type of computation 

(E- or H-polarization,backscatter or bistatic echo width), generate the desired data and 

f 

store the resulting output. Some initial post-processing is also performed. For instance, 
if the near-field values on the grid are stored (this option is only available for bistatic 
computation), an additional file may be generated for a contour plot. 

The following menu is produced by FECGFFT during its execution. 

****** Main Menu ****** 

Pre-processing 

(1) Load Finite Element Mesh File 

(2) Set new CG parameters 
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Analysis 

(3) E-polarization (Backscatter) 

(4) H-polarization (Backscatter) 

(5) E-polarization (Bistatic) 

(6) H-polarization (Bistatic) 


Post -processing 

(7) Generate 3-D plot file 


Test Routines 

(10) Element node ordering 

(11) Test integral matrix: scattered fields 

(12) Free-space field comparision 

(20) quit 


Item (1) allows the user to load a mesh data file generated from one the mesh gen- 
erators to be discussed later. Actually, any mesh generator may be used, but the file 
must contain the correct information and format. This information can be found by 
examining the module MLOAD. 

Item (2) allows the user to change the CG residual error tolerance value and interval 
for printing the iteration number and the associated residual error. 

Items (3) and (4) are selected for the generation of backscatter data for E- and H- 
polanzation, respectively. When either of these is selected, the starting angle, stopping 
angle and angle increment will be prompted. The file name for the far-field data is also 
requested. A response of “none” will produce no file. A prompt for the pad size will then 
be requested. The suggested response of “1” will automatically determine the proper 

pad size. When the program has finished generating the desired data, it returns to the 
main menu. 

Items (5) and (6) are selected when bistatic data for E- or H-polarization are desired. 
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When either of these is selected, the incident angle will be prompted followed by the 
starting angle, stopping angle and angle increment. The file name for the far-field data 
is requested, followed by a file name for storing the nodal field values. A prompt for the 
pad size follows as before. 

Item (7) allows the user to generate data in MPLOT format for contour plots. At 
the present time, only rectangular bodies will work for this option. 

Items (10)-(12) direct the user to test routines, not used for normal operation. 

The pertinent files which contain groups of subroutines associated with the accom- 
panying description are as follows: 


| file name 

brief description 

fe.cgfft.ftn 

main program 

fe_vec_sub.ftn 

vector operation subroutines 

feJo.sub.ftn 

file i/o routines 

fe_test3_sub.ftn 

near-field test routines 

fe.cross^sub.ftn 

cross-term subroutines 

fe_matrix_sub.ftn 

FE matrix routines 

fe_test5_sub.ftn 

node order test routine 

feJft_sub.ftn 

FFT routines 

fe_three-sub.ftn 

three-dimensional plot data generation 

fe_field_sub.ftn 

near /far field computation 

fe_ftest.sub.ftn 

free-space test routine 

fe_summary_sub.ftn 

generates session summary 
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These files must be compiled and linked prior to execution. 

7.2 Mesh Generator for Curved Bodies 

This program is under development for various specific types of bodies. It may, 
however, be used to generate a mesh for virtually any desired body. The mesh generation 
is accomplished by first dividing the region between the impenetrable surface (if any) and 
the rectangular enclosure into first, second or third order serendipity elements. These are 
subsequently mapped to a square domain, subdivided and mapped back. Examples of 
these are shown in fig. 7.1. It works well for modelling regions with curved boundaries, 
but generally produces more unknowns than necessary for the solution method. 

An input file to this program may be generated either with option (2) from the main 
menu, or manually. Selection of option (3) from the main menu processes this file and 

places the results in a specified output file. The output file is then used as input to the 
program FECGFFT. 

Running the program produces the following menu: 

Main Menu 

(1) Preferences 

(2) Create an input file 

(3) Process an existing input file 

(5) Plot routine 
(10) Quit 

Item (1) has not been incorporated as of yet. The selection of Item (2) produces the 
menu: 
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Figure 7.1: Typical serendipity elements used in the region descretization process 
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Input File Creation Menu 

(1) Conducting Strip 

(2) Rectangular Coated Slabs 

(3) Coated Ogive 

(4) Circular Cylinder 

(10) Return to Main Menu 


Only Items (3) and (10) are operational at this time. Selection of (3) yields the menu: 


Ogive Menu 


(1) Enter geometry 

(2) Modify geometry 

(3) RETURN 


Choosing Item (1) results in a series of prompts outlined as follows: 

1. a,b for the coating, where a = height of the arc and b = half-length as indicated 
in fig. 7.2 

2. the relative permittivity and permeability of the coating 

3. a,b for the conductor 

4. sampling interval (in wavelenths) at the integration boundary 

5. number of circumferential samples in the free-space region 

6. number of circumferential samples in the material 

7. arc distance along coating 

8. arc distance along conductor 

9. comment to appear in the input file 
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Figure 7.2: Arc used for part of an ogival structure. 

10. input file name 

Upon completion of the session, a file is generated to be used at the input to FECGFFT. 
Selecting Item (3) results in a prompt for the input and output filenames. 

Manual generation of the input file requires that the scattering body be surrounded 
by a rectangular boundary displaced approximately one element from the body. The 
region within the rectangular boundary and the impenetrable body surface (if present) 
is then subdivided into either linear, quadratic or cubic elements, examples of which are 
given in fig. 7.1. Note that every node and side of each element is numbered as indicated. 

The output file contents are listed as indicated in tables 7.1 and 7.2 with the variable 
and descriptions in table 7.3. The first four lines are self explanatory. The next group 
of lines contains the coordinates of the nodes, and the order of these pairs determines 
the global node numbering scheme. Two real numbers followed by a “C” are assumed 
to be in cylindrical (r,0) coordinates centered at the previously specified value. An “N” 
following the “C” will change the center coordinates to x c ,y c . Immediately after the 
node coordinates definition, the elements or “blocks” are defined. The local/global node 
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relationship defines the block. This format is repeated for each of the N e elements. 
The elements must then be connected by specifying the sides of adjacent elements. This 
avoids the time-consuming task of comparing the coordinates of every node to the others 
for spatial commonality. The impenetrable and integration boundaries designation are 
present for a similar reason. Finally, the material properties are requested. The order of 
these determines the number to be used in determining the element material property 
number. The free-space value is always present. 
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variable 

type 

description 

N e 

I 

total number of modelling elements 

fi 

I 

put element numbers on mesh (l=yes, 2=no) 

f 2 

I 

put node numbers on mesh (l=yes, 2=no) 

h 

I 

put material numbers on mesh flag (l=yes, 2=no) 

u 

I 

surround mesh with scale (l=yes, 2=no) 

fs 

I 

generate PostScript file (l=yes, 2=no) 

N n 

I 

total number of nodes 

(x,y) 

R,R 

cartesian coordinate pair 

( r >9) 

R,R 

cylindrical coordinate pair 

OcS/c) 

R,R 

coordinates of arc center 

Mi 

I 

material number of t’th element 

N < 

I 

number of samples on side i of element j 

Oi 

I 

order of ith element 

L ~ G(i) 

I 

global node number of ith local node 

N c 

I 

total number of element sides in contact 

e, 

I 

element number 

3 i 

I 

side i of element j 

N bc 

I 

number of elements adjacent to conducting boundary 

N bi 

I 

number of elements adjacent to integration boundary 

N p 

I 

number of entries in material table 


Table 7.3: Variable description and FORTRAN declaration type. 
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Figure 7.3: The mesh of a conducting circular cylinder. 


Example 

The following file is an example of a perfectly conducting cylinder of radius .5 A in 

free space. A figure of the resulting mesh is shown in fig. 7.3. 

Circular cylinder 8 Aug 1989 

8 , 1 , 0 , 

2 , 2 , 2 , 2, 2 

40, 

•5,90. ,CN,0. ,0. 

•5,112.5,C 
.5,135. ,C 
•5,157.5,C 
.5,180. ,C 
. 5,202 . S ,C 
.5,225. ,C 
. 5 , 247 . 5 , C 
.5,270. ,C 
-5,292.5,0 
• 5,315. ,C 
.5,337.5,0 
.5,0. ,C 
.5,22.5,0 
.5,45. ,C 
-5,67.5,0 
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.525,90. ,C 
.64,135. ,C 
.525,180. ,C 
.64,225. ,C 
.525,270. ,C 
.64,315. ,C 
.525,0. ,C 
.64,45. ,C 
.55, .55 
.275, .55 
0. , .55 
-.275, .55 
-.55, .55 
-.55, .275 
-.55,0. 

- . 55 , - . 275 
-.55, -.55 
-.275, -.55 
0. ,-.55 
.275, -.55 
.55, -.55 
.55, -.275 
.55,0. 

.55, .275 

1- ST BLOCK 

0 , 0 , 1 , 

H. 3,8, 

I . . 1 ., 

I, 15,25,27,16,24,26,17 

2- ND BLOCK 

0 , 0 , 1 , 

II, 3,8, 

1-,1., 

3,1,27,29,2,17,28,18 

3- rd BLOCK 

0 , 0 , 1 , 

3,11,8, 

1..1., 

31,5,3,29,19,4,18,30 

4th BLOCK 

0 , 0 , 1 , 

3.11.8, 
l-.l., 

33,7,5,31,20,6,19,32 
5th BLOCK 
0 , 0 , 1 , 

11.3.8, 

I. ,l., 

33,35,9,7,34,21,8,20 
6th BLOCK 
0 , 0 , 1 , 

II, 3,8, 
l.,l., 

35,37,11,9,36,22,10,21 
7th BLOCK 
0 , 0 , 1 , 

3,11,8, 
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1 , 1 ., 

11 , 37 , 39 , 13 , 22 , 38 , 23,12 

8th BLOCK 

0 , 0 , 1 , 

3 , 11 , 8 , 

1 . , 1 . , 

13 , 39 , 25 , 15 , 23 , 40 , 24, 14 
CONNECTION 
8 

1 . 1 . 2. 3 

2 . 1 . 3. 4 

3 . 2 . 4. 4 

5 . 1 . 4.2 

6 . 1 . 5. 3 

7 . 2 . 6. 3 

7 . 4 . 8. 2 

8 . 4 . 1.3 

Perfectly Conducting Boundary 

8 , 

1.2 

2,2 

3.3 

4 . 3 , 

5 . 4 , 

6.4 

7.1 

8.1 

Integration Boundary 

8 . 

1 . 4 , 

2.4 

3.1 

4.1 

5.2 

6.2 

7.3 

8.3 

Dielectric property table Epsr, Epsi, Mur, Mui 

1 . , 0 . . 1 . , 0 . 


The program contains the following files: 
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j file name 

brief description |j 

gen .ft n 

main program 

gen_sub.ftn 

contains associated subroutines 

fe_grid_ttz_sub.ftn 

for plotting a mesh with rectangular 
elements on the Apollo screen using 
graphics primatives 

fe_grid_sub.ftn 

for plotting a mesh with triangular 
elements on the Apollo screen using 
graphics primatives 

fe_post_sub.ftn 

for generating a postscript version 
fe_grid_sub.ftn 


These programs should be compiled with the SAVE option and linked before execution. 

7.3 Mesh Generator for Rectangular Bodies 

This program is useful for generating the mesh associated with coated rectangular 
bodies. Executing the program produces the following menu: 

**** Mesh Generation Menu **** 

(1) Conducting Strip 

(2) Composite Bodies 

(3) View an existing file 

(10) Quit 

Item (1) should be ignored, since the data file it generates for the strip does not 
distinguish between the nodes above and below the strip. It is thus invalid for H- 
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polarization computations. Item (2) allows the user to generate a rectangular composite 
body. Upon its selection, the user is promped for the following items: 

1. the size of the square building block cell (in wavelengths) 

2. the permittivity and permeabilities of the various material layers for the structure 

3. the length and width of the main scattering structure in integer multiples of the 
initially specified building block size in 1. 

4. the type of scattering body (conductor or material) 

5. the number of layers for each side plus the material number from the material table 
generated in 2. 

6. the number of cells between the scattering body and each of the four boundaries 
(usually 0, unless a larger grid is desired) 

7. the name of the output file to be used by the FECGFFT program 
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Figure 7.4: The mesh of a rectangular partially coated cylinder. 


Example 


To generate a .5A x 2A conducting body with a .1A material coating of e r = 5. - j. 5 

and = 1.5 - 3 . 1 ove, aide. (1) and (2) (sea fig. 7.4, the following ootpu. is displayed: 

**** Mesh Generation Menu **** 

(1) Conducting Strip 

(2) Composite Bodies 

(3) View an existing file 

(10) Quit 

2 

Enter del (size of building block) in wavelengths 


Enter dielectric materials to be used [(-1..0.) to quit) 
(Remember, Imaginary parts <*0.) " ^ 


( 1 . 000000 , 0 . 0000000) 
( 1 . 000000 , 0 . 0000000 ) 


Epsilon 1 * 

Mu 1 

Enter Epsilon 
(5. ,-.5) 

Enter Mu 2 
(1.5,-. 1) 

Enter Epsilon 3 

(- 1 .. 0 .) 

4 Enter length and width of main body in units of del 
Q Main body composition: (0) Conductor (1) Dielectric 



1 


I 

2 1 |4 

l___ I 

3 

Enter Onumber of dielectric layers sides 1, 2, 3, 4: 

2 , 2 , 0 ,0 

Index Epsilon Mu 

1 1.000 0.000 1.000 0.000 

2 5.000 -0.500 1.500 -0.100 

^Material property number for: Side 1 layer 1 

^Material property number for: Side 1 layer 2 


Index Epsilon 

1 1 . 000 0 . 000 

2 5.000 -0.500 


Mu 

1 . 000 0 . 000 
1.500 -0.100 


^Material property number for: Side 2 layer 1 

^Material property number for: Side 2 layer 2 

^ Number of rows and columns of blank cells surrounding the body 

^Generate PostScript file? (l*yes, 2®no) 

Enter file name for data storage 
test_out ° 


**** Mesh Generation Menu **** 

(1) Conducting Strip 

(2) Composite Bodies 

(3) View an existing file 

(10) Quit 


The output file test.out contains the mesh information required by FECGFFT. 

Item (3) provides for viewing the plot on an Apollo screen. Upon its selection, the 

user will be prompted for a file name. Entering the interactive mode results in the 
following menu: 

Interactive Mode Menu 
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(1) Max and Mins 


(2) Picture orientation 

(3) Picture size 

(4) Picture offset 

(5) Tick spacing 

(6) Legend contents 

(7) Legend offset 

(8) Legend label size 

(9) Label contents 

(10) Label size 

(11) Number size 

(12) Number format 

(13) Print option flags 

(14) View on screen 

(15) Get hard copy 

(16) PsPreview hard copy 

(17) Reset default values 
(20) Return to main menu 


Currently, options (6)-(10) have not yet been incorporated. The remaining items 
self-explanatory. 

The programs contains the following files: 


file name 

brief description 

mgen_lin_nc_new2.ftn 

main program 

fe_grid_sub.ftn 

for plotting a mesh with triangular 
elements on the Apollo screen using 
graphics primatives 

fe_post.sub.ftn 

for generating a postscript version 
fe_grid_sub.ftn 


which should be compiled with the SAVE option and linked before execution 


are 
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